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^ Abstract 

Based on the analysis of a certain class of linear operators on a Banach space, we provide a closed 
form expression for the solutions of certain linear partial differential equations with non-autonomous 
input, time delays and stochastic terms, which takes the form of an infinite series expansion. 
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S Introduction 

Linear differential systems are ubiquitous in pure and applied mathematics, either as models, approxima- 
CN tions, but also because the stability of solutions of nonlinear differential systems reduces to the study of 

^ linear systems. Such systems might include stochastic terms (see [31), temporal delays (see HI), and also 

ly-j encompass the case of partial differential equations. Apart from the simplest linear finite-dimensional differ- 

OO ential equations, finding closed forms expressions for the solutions of general linear differential systems is 

very complex. In this paper, based on the treatment of evolution equations as algebraic equations in a suitable 



^ Banach space, we propose a closed-form expression for the solution of linear, non-autonomous, stochastic, 

time-delayed partial differential systems. Application of this framework to several classical examples such 



X 



as the delayed Ornstein-Uhlenbeck process or the stochastic heat equation are developed in sections 3.3 and 



3.4 This expression is especially useful to understand the dynamics of weakly connected linear learning 



neural networks, problem which motivated the development of this more general framework (sections [T] and 



3.51. 



1 Motivation: neural networks with learning 

Consider a neural network made of n G N* neurons described by their membrane potential at time t written 
V{t) S W^. The connectivity matrix W S W^^"' contains the weights of the connections between each pair 
of neurons. Learning in biological networks corresponds to the slow evolution of the connectivity. Neural 
network with Hebbian learning can be modeled by the very simple coupled model: 

eV =-lV + W-V + Iit) 

W = -kW + V(^V ^ ^ 
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where the positive parameters e governs the ratio between the potential and the connectivity time scales, / 
the leak conductance and k the learning ratio, and I{t) G is a r-periodic external input. The symbol 
(g) denotes the tensor product: {V ® V^ij = A companion paper (in preparation) proves the well- 

posedness, and analyzes the dynamics and the equilibrium connectivity of such systems in the limit e — )• 0, 
i.e. when learning occurs on a much slower scale than the activity. In that asymptotic regime, the connectivity 
is well approximated by the equation W = —kW + ^ jQ^wi^) ^ V^^r{t)dt where V^r is the periodic 
solution ofV = -IV + W - V + I{t) where W is assumed constant. Seeing V^i G C^([0, t[, M") as a 
"semi-continuous" matrix of size n x [0, r[, such that V^r^^ = (Vjy leads to the formulation 

T 

The characterization of the solutions of these equations is therefore important to understand what is 
learned by such networks. The following results will provide a way to compute the equilibria as shown in 
section [331 

2 Framework and General Result 

The framework we develop here is based on extending notions of matrix calculus to infinite dimensional 
spaces. The linearity of the equation motivates to extend some finite-dimensional linear algebra and matrix 
concepts to infinite-dimensional spaces. 

We consider in the manuscript linear equations in a Banach space C of real functions of time t and a 
variable x £ E, called space variable, where E can either be a finite set {!,..., N} (in which case C is 
equivalent to the space of M^-valued functions), countable or continuous, typically M, in which case C is 
a space of two-variables functions. The particular problem under consideration governs the choice of the 
space C, in particular including regularity or integrability properties (typically C is a or a Sobolev space). 
Similarly to a matrix notation, we denote the value of X £ C at {x,t) £ E x Rby X^t- 

Let E denote the space of bounded linear operators on C. We are interested in solving equations of type 
CX = B where C ^ £ (this operator may involve differentials in time and/or space) and B £ C. We 
will restrict the study to a class of operators of a particular form we now detail. To this end, we introduce 
two kinds of linear operators on C: the space operators L acting on the first (space) variable, i.e. linear 
operators on M^. If E is finite, this set is reduced to the matrices. If E is equal to M.'^, it contains all the 
linear operators acting on functions of the space variable, in particular, under suitable regularity conditions, 
integral or differential operators. The action of the space operators L on a function X £ C is denoted L • X 
(acting on the left). The time operators essentially act on the second (time) variable, and the transform 
might depend on the space variable x. In other words, these transforms TZ can be represented by a family of 
operators {Tlx ,xGE) such that for any x, Tlx is a linear operator on (M) . The action of a time operator Tl 
on X G C is written X ■ Tl (acting on the right). In the paper, we will mainly be interested in diagonalizable 
time operators. Diagonal operators in the time domain are operators Tl whose action can be written in the 
form {X ■ TVjxt = f{x, t)Xx,t- This class includes for instance all linear differential time operators, which 
are diagonalizable in the Fourier basis. Another class of time operators we will be considering is the class 
C of convolution operators with respect to time. Given a finite measure g of M, the convolution operator 
Tg ^ C associated with g is defined as [X ■ Tg)^^ = Xx(^t~s)dg{s). Such operators are generalizations 
of Toeplitz matrices generated by g, with, loosely speaking, infinitely many rows and columns. An important 
property of the convolution operators is that they are diagonal in the Fourier basis. 
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For L a space operator and TZ a time operator, we define tlie Kronecker product L (g) 7?. as tiie mixed 
operator of £ sucii that (L ® TZ) {X) = L ■ {X - TZ). Note tiiat tiie product becomes associative wiien 7^ is a 
convolution operator which will be the case in section[3] This definition extends the property of vectorization 
of the Kronecker product of matrices in linear algebra (see e.g. fTl). 

The main technical result of the paper is given in the following: 

Property 2.1. Let C = A B + Idc ® T) be a linear operator, for A a space operator and B, T) co- 
diagonalizable time operators, with B invertible. For the sake of simplicity, we assume that they are diagonal 
in the natural time basis, and denote for x G E, Bx = diagti=M.(t){x,'t)^ andT>x = diagt^^{d{x,t)^. We 
assume that va.ix,t \ b{x, t)| > and the spectral condition: 

def Www 

3 / € M* such that A = f — , , ^ < !> (2) 



where W = IM^e + A and W^W = supxj^o ^^^x\\^^ operator norm. Then A B + Idc ® ^ 

invertible and its inverse reads: 

A^B + Idc0Vy =-^W''^d^ag,,J .(.,),,^J - (3) 

Remark The spectral condition is merely a technical sufficient condition for the convergence of the series. 
The relatively formal setting and assumptions will become clearer in the applications, section^ 

Proof. The direct introduction of the inverse can appear artificial at first sight. However, this formula is a 
natural extension of the discrete-time case where direct linear algebra and Kronecker products calculations 
quite simply provide a closely related expression the interested reader can readily derive. 

In order to prove the proposition, we first need to prove that the operator indeed exists, and that it 
constitutes the inverse of C. It is easy to show that under the assumption of the proposition that the sequence 

of operators in £ defined by: Mn ^= —Ylk=o^^ ® ^^'^^^'^^{, b{ • t){i constitutes a Cauchy 

sequence in £. Since C is a Banach space, so is £, and hence the sequence {Mn)n converges. The limit of 
this sequence is our inverse candidate, and is denoted as the infinite series Q. 

In order to prove that this limit is indeed the inverse of C, we compute the limit of (Mj^ ° C)X (or 
similarly (£ o Mn)X) for a given X £ C.lt is easy to show, developing the series, that we have: 

N 

o c)x)xt = -y-w'. -^-^ + • , . = Xxt - w"^^' ^ 



b{.,t)j V ''(•'*); V ''(•'*); 

where for y G C denotes the application E such that y,t(x) = Yxt- Here again, the assumptions 
of the proposition ensures that the second term vanishes as goes to infinity. □ 



3 Application to solving linear time-delayed Stochastic Partial Differential 
Equations 

In this section we make explicit the use of the inversion formula ([3]) in the case of linear delayed, stochastic, 
partial differential equations. Several examples with different convolution operators will illustrate the main 
result of the section stated in theorem [XT] 
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3.1 General Result 

Let be a Hilbert space, typically M" for n S N, L^(M") or a Sobolev space of applications on M". We 
consider a probability space {Q, T, P) satisfying the usual conditions and B a standard adapted ^Y-Brownian 
motion (for the existence and properties of this object in infinite-dimensional spaces, see |i4j Chapter 4]). We 
aim at solving the non-autonomous time-delayed stochastic differential equation: 

dX = {A - {X * g) + I) dt + T,- dB 

= Co ^ LliR*,) (4) 

with : X ^ X linear, / G C(]R-|-, A') an external input, g a finite measure of the real line supported 
on ]R_|_, i.e. a causal measure, and * denoting the convolution. Existence and uniqueness of weak solutions 
for such equations is ensured, see e.g. jBllH. We consider the case where the system has a unique strong 
solution. In the case where X = M", this occurs under the assumptions of the section, see e.g. |3, Chapter 
5], and in the infinite-dimensional case, we need to assume that i? is a genuine Wiener process (i.e. the trace 
of the covariance matrix is finite, and the initial condition is in the domain of A, see [4|). The solution of this 
stochastic differential equation at time t e M+ is defined by the integral equation X{t) = Co(0) + /q [A ■ 
{X * g){s) + /(s)) ds + S • dB and = Co- 

This problem can be set in the framework described in section [2] using a transformation inspired by the 
classical Fourier transform of the solution in the time domain. To perform this transformation rigorously 
in our particular stochastic setting, we stop our processes at a finite time r > 0. We define X^ : t G 
M — l[o the restriction of X to the compact support [0,r] and null elsewhere. Similarly, define 

Ir = 1[o.t]-^ ^rid dBr = l[o We have: 

Theorem 3.1. For all r G M_|_, choose I G M* and W a space operator such that W = lldc + A. If 
the spectral condition Q is satisfied, i.e. in the present case \\W\\ < infg{|Z + then the solution of 

equation (|?]) is given by 

+00 

Xr = ^W''- (Co(0)5o + /r + S • dBr^ ■ U ■ (5) 
fc=0 

where U = F- rf^o^^eR • -^"^ V = F- diag^eM.(^j^(0!^') ■ and 7^ = + A • (Co * g)- 

The notation [dBr - U) ^ stands for the stochastic integral JJ"*"(*'^) UstdB{s) which is square integrable on 
[0,r[. 

Remark: r/te convergence of the series ([5]) occurs as soon as the spectral condition ([2]) is satisfied on the 
subspace spanned by W'^ ■ (Co(0)(5o + /r + S ■ dBr) -U -V^. 

Proof. First, note that A ■ {X * g) = A ■ {Xj- * g) + A ■ {Qq * g) yielding the equation on Xr'. dXr = 
[A • {Xj- * g) + It) dt + T, ■ dBr. Thus, the initial condition on X acts as an external input on Xr. In the 
deterministic finite-dimensional case, it is well known that differential operators are diagonal in the Fourier 
basis. Based on this result, we introduce the Fourier transform F of equation (|4]l for a fixed (x) G il. As 
mentioned, for almost all G il, the processes involved are bounded, hence the function of time, on the 
compact interval [0, r], is square integrable in time. Let : t G M — )• e^"^'^"^^^ Xr{t) for C G M the Fourier 
variable and X is the unique solution of equation [4] Ito formula yields for t < r 

dZ^{t) = - 2iTiiZ^{t) + A • (Z« * g){t) + e-^'^*^Ir{t)jdt + e'^^^^^S • dBr{t) (6) 
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Let us denote by : ^ € M £ Z^{s)ds the Fourier transform of Xr and U : i ^ /J" e-^^'^*^ir{t)dt. 
The process Br is the well-defined stochastic integral /J^ e^'^^'^^^dB{t). The integral form of equation ([6]), 
using the fact that the convolution is diagonal in the Fourier basis, denoting D = diag^gK ( — 2i7r,^) and 
Q = diag^^]R[g{^)), leads to the functional equation: 

Z-{t)- Z-{0)= A-Xr-G + Xr-D + ir + ^-Br 



Applying proposition 2.1 for a fixed u ^ 0, where C is the set of square integrable functions on [0, t[ which 




is a Banach space, we obtain: 

+00 ^ 

Xr = ^W''.[-Z-{r) + Z-(0) + ir + ^-Br)- d^ag^,u{ ) (7) 

k=o 9[c,){i' + ) 

We now take the inverse Fourier transform of this expression by applying the time operator T^^. First of all, 
we perform the inversion on the terms 1^ = Ir ■ T . It is easy to show that 1^ ■ diag[—— — liTv^^k+i ) " ~ 

It - U ■ V'^. Similarly, the term (^Br ■ diag[ ^^^^^^^1^^^^^^ ) ■ T^^^ can be written dBr - hi -V^. 

Moreover, for x £ {0, r} an easy computation shows that (^Z'{x) ■ diag(^ -(^^^(^i^\'^i^k+i ) " = 

(^X{x)6x) • U • V^. Furthermore, the operators U and V are causal, i.e. if Y has a support C [c, +oo[ then 
Y -U ■ V'^ also has a support C [c, +oo [. Indeed, n : ^ i— ;g(g)j^2t7rg corresponds to the transfer function of 
a closed loop filter shown on the right, and hence U is clearly causal since g is. 

V is also causal as the convolution of g and U. This implies that the contri- 
bution of Z{t) vanishes in equation (j7]l since it has its support in [r, oo]. 

3.2 Computational Remarks 

Truncations of the formula (|5]) provides approximations of the solution of system Q. We observe that the 
smaller the difference between the spatial operator A and a multiple of the identity, the more accurate a 
truncation of this expansion. In other words, this expansion is particularly useful if A is a small perturbation 
of the (scaled) identity (e.g. the case of weakly connected linear neural networks). 

This representation allows development of new numerical schemes for the simulations of the solutions 
of system Q. For simplicity, consider the case where E = {!,••• ,n}. To approximate the solution over 
the interval [0, r] define a time step At an number of points T = r/At G N and replace U and V by the 
Toeplitz square matrices U and V, generated by f € {0, • • • , T — 1} i— )• j/^^^^'^* u{s)ds, where u is the 
function generating U (and similarly for V). The number of operations needed is 0{{k + l)nT{n + InT)) 
since the product with a Toeplitz matrix, as a convolution, has a cost O(TlnT). This scheme has a first 
order accuracy, 0(dt^ + (ix + A*^+i) where 7 is equal to 1 for deterministic equations or ^ if stochastic. In 
comparison, the Euler-Maruyama method has a complexity of 0{T{n^ + n where 9 is the support 

of g and an accuracy of O [df^ + dx) , comparable to the expansion method in both aspects. 

Two interesting advantages of the expansion over Euler-like methods are that (i) it is parallelizable and 
(ii) it appears to be numerically very stable, i.e. large At do not lead to a diverging scheme. 

3.3 Examples 

Let us now treat some classical problems that are solved in the present framework. 
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• Ornstein-Uhlenbeck process: The simplest example is the Ornstein-Uhlenbeck process with no de- 
lays (i.e. g = 5q). In that case, ^ = 1, and therefore, for all / G M, inf^{|/ + 2i7r^|} = |/| and 
the expansion is valid if there exists / G M* such that ||/ + ^|| < |/|, i.e. for any operator A whose 
spectrum is bounded and entirely contained in the left or right half plane. For negative matrices A 
(i.e. / >0)7h=Z^ = Visa Toeplitz operator generated by the function h : z ^ e'^^^H{z) 
with H the Heaviside function. Therefore, the solution oiX = A- X + I + T, - dB can be written as 

= Ylt=QW^ ■{Xo5o+lT+'^-dBr)-T^"^^ .If A Idc {t.g. in one dimension), the terms for /c > 
vanish in the above equation and we get the simple well-known expression Xr = {Xo6o+lT+dBr)*h. 

• Exponentionaly distributed delays: Let us now treat the case g : x ^ [3e~^^ H{x). In this case, 
g{^) = ^^2CTg • Therefore, |^ = — 27r(^^^ — z^) which corresponds to the red curve in the left 
picture of figure [T] Operators A satisfying the spectral condition |2] are the ones whose spectrum is 
contained in an open ball centered at — / that does not intersect the red curve (blue disks of in figure 

When A is negative, the operators U and V can be made completely explicit. Indeed, observing that: 

= (g^,,,g).„g^ = ip--^+2^-.ml^+2^.i) A = ^1^4^, the Operator V is the 
convolution operators generated by j3{h^*hj^^ with /i-t : 1 1— )■ e~^^2~^H{t). Similarly is generated 
by /3(/i„ * hj^ + ^h'_ * hj^) . Even more explicitly, for /3 > 4/, V is generated by i i— ^e~2^sh[^t) 

and Ubyt>-^ se"^*(| sh{^t) + /S.ch[^t)^. When 4/ > /?, a similar result holds replacing the 
hyperbolic functions ch and sh hy cos and sin. 

• Single fixed delay: For g = 5q + ade, we have g{^) = 1 + ae^*'^^^. The convergence domain of 
the expansion (condition [2]) is shown in the middle and right pictures of figure [TJ for two different 
a G M+. The red curve seems to be living on the 2-dimensional projection of a simple 3-dimensional 
cone of revolution whose section is a circle. In that case, it appears quite difficult to express U and V 
in a simple form, though their Fourier transform is explicit. 

Remark: As illustrated in the previous example, a procedure to find the constant / such that the expansion 
converges consists in plotting on the same figure the complex eigenvalues of the spatial operator and the red 
curve ^ G M I— )• |^ G C related to the time operators. If there exists a ball centered on the real axis which 
contains all the eigenvalues and that does not intersect the red curve, then choosing — / as the value of its 
center ensures that the expansion will converge. 

3.4 Stochastic heat equation 

Let us now deal with a classical the stochastic heat equation on (described as the interval [0, 1] where 
and 1 are identified) as a classical example of linear partial differential equations: 

du 

— (x, t) = An(x, t) + 11(2;, t) + cr?7(x, t) 

with periodic boundary conditions, where A is the Laplacian on S^, v is an external forcing and r/ is a 
multidimensional white noise. The input v{x, t) is set to 5x=xo{x) and we take the initial condition u{., t = 
0) = 0. 

The Laplacian operator has eigenvalues —An'^k'^ with G N, corresponding to the eigenvectors cos(2fe7rx) 
and sin(2A;7rx). Since the eigenvalues are not bounded, it is not possible to find a suitable constant I to define 
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Figure 1: The three different pictures correspond to different time-convolution kernels g. The red curves 
are the parametric plots of ^ S M i— )• G C and the blue balls are examples of sets within which the 
eigenvalues of the space operator A need to live for the expansion to be well-defined. The eigenvalues have 
to be contained in a single ball. The center of each ball is — /, for different / G M. To satisfy the spectral 
condition (|2]) the balls cannot intersect the red lines, (left) Exponentialy distributed delays with /3 = 27r and 
i G [-20, 20]. (middle) Single delay with a = 2 and 6* = 1 and ^ G [-20, 20] (right) Single delay with 
a = 0.3, = 1 and^ G [-5,5] 



the solution of the heat equation in our framework. However, the semi-discretized in space equation over- 
comes this problem by preventing the existence of very fast oscillations (corresponding to large eigenvalues 
of the Laplacian). We choose to discretize the space with N points regularly spaced, corresponding to a 
discretization step dx = The resulting equation corresponds to Q in dimension N, with g = 5{t) 

and A G M^^^ such that An = -2/dx^, Aij = Xjdx^ if i = j ± 1, A^n = Ani = l/dx^ accounting for 
the periodicity of the medium, and Aij = otherwise . This matrix has eigenvalues in [— 4/(ix^,0]. This 
suggests the choice / = l/dx^ so that all the eigenvalues are in this ball a center — / and radius I. This ball 
intersect the imaginary axis only in (corresponding to spatially constant functions), so convergence issues 
might arise if one of the terms W'' •(/,- + S • dBr) • ^ • is spatially constant, which clearly never occurs in 
our case. Therefore, our expansion is well-posed and provides a numerical scheme to compute the solution 
as shown in figure [2]d. In that figure, we exhibit the fact that the solution is well retrieved by the expansion, 
and the error compared to Euler's scheme with a time step dt = 0.01 (Fig. [2]b) is more than two order of 
magnitude smaller than the solution. An interesting point of this method is that it works for any time step 
interval dt which is not the case for the Euler method which rapidly diverges as soon as the CFL condition 
is not satisfied for instance. Moreover, extending the approach to a delayed formalism g ^ 5 i?, costless in 
our framework. 



3.5 Weakly connected networks 

Let us now go back to the problem of weakly connected neural networks introduced [T] The application of 



theorem 3.1 to the Ornstein-Ulhenbeck process (see first example in section 3.3), allows characterizing the 
equilibrium connectivity of system ([T]) with a null initial condition as the solution of: 



1 - M-^ 
K Ik+i ■ 

k,q=0 



+1' 



. w- 



iq+l 



In the weakly connected limit, i.e. 



+00, it can be shown that \\W*\\ < I and we can find a good 
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Figure 2: Application of the expansion method to the stochastic heat equation on the circle with a Dirac 
source on the neuron in the middle, a) Space-time diagram of the solution given by the expansion method 
for dt = 0.01. b) Space-time diagram of the error between the solution in a) and the solution given by Euler's 
method, c) Space-time diagram of the solution given by the expansion method for dt = 1. d) L2 norm of 
the terms in the expansion. The parameters for these simulations are n = 100, number of time steps = 500, 
cr = 0.1 and / = 2. 



approximation of W* in the form 

^* = §^-r..'7;'./' + o(§) 



where 1^ = sup^gjj II2. This expression makes explicit the fact that the equilibrium connectivity stores 
not only the spatial, but also the temporal correlation of the inputs. The method also makes it possible to 
go at arbitrary order and consider also correlated forcing noise together with delays in the communication 
term. 
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